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Modeling of cotunneling in quantum dot systems 
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Transport through nanosystems is treated within the second order von Neumann approach. This 
approach bridges the gap between rate equations which neglect level broadening and cotunneling, 
and the transmission formalism, which is essentially based on the single-particle picture thereby 
treating many-particle interactions on an approximate level. Here we provide an alternative presen- 
tation of the method in order to clarify the underlying structure. Furthermore we apply it to the 
problem of cotunneling. It is shown that both elastic and inelastic cotunneling can be described 
quantitatively, while the transmission approach with a mean-field treatment of the interaction pro- 
vides an artificial bistability. 

PACS numbers: 73.23.Hk,73.63.-b 
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I. INTRODUCTION 

Nanostructure technology allows for the fabrication of 
small structures, such as quantum dots [1], nanowires 0, 
carbon nanotubes Q, as well as molecular structures, 
whose electronic properties are dominated by a small 
number of electrons. The insertion into electrical cir- 
cuits exhibits current-bias {IV) relations, which strongly 
differ from conventional resistors. Typical features are 
the suppression of current for low bias due to Coulomb 
blockade and pronounced current steps, whose position 
can be easily controlled by a gate bias, thus suggesting 
transistor action. 

Due to the small spatial dimensions, the level quantisa- 
tion with a spacing /S.E is of fundamental relevance. Fur- 
ther relevant energy scales are: The Coulomb interaction 
between electrons (or similarly any other type of many- 
particle interaction [3]) of the order U = e^/C, where 
C is the geometrical capacitance; the tunneling rate T /h 
for the transition of particles between the structure and 
the contacts; as well as the thermal energy kBT. For 
typical nanostructured electronic systems studied exper- 
imentally, these energies are all in the range of 0.1-10 
meV, where F can be even smaller. Generic examples of 
/^-characteristics are shown in Fig. [T] 

The task to describe the current through such struc- 
tures quantitatively is an evolved topic and a large variety 
of approaches has been used within the last two decades. 
The general starting point is to divide the Hamiltonian 
of the system as 
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Here Hd describes the central region (such as a single 
or multiple quantum dot structure, a molecule, a car- 
bon nanotube, or a nanowire), which we refer to as the 
quantum dot in the following. 
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FIG. 1: Schematic current-bias relation for a quantum dot, 
showing the qualitative influence of the key energy scales. 



describes the leads, where a =tji denotes the spin, k 
labels the spatial wave functions of the contact states, 
and £ denotes the lead (typically i = L/R for the left 
and right lead, respectively). In most works a thermal 
occupation of the leads 
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with chemical potential ii£ is assumed, where the bias 
eVbias = tJ-L — fJ-R drives the current in a two termi- 
nal system. Conventionally noninteracting leads are 
assumed, [s?] Finally, Ht describes the coupling between 
the leads and the quantum dot. 

If Hu does not contain many-particle interactions, it is 
straightforward to derive a transmission formula (see @ 
and references therein) yielding a particle current from 
lead £ into the quantum dot 
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j£ = 2^ y Ti^i,{E) [ME) - M{E)] 
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Here the transmission function is determined by the 
single-particle spectrum of Hu in connection with the 
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coupling Ht, which essentiahy broaden the levels. A for- 
mal tool is the nonequilibrium Green function formal- 
ism An inclusion of many-particle effects within 
the mean-field level is straightforward and exchange- 
correlation effects can be taken into account by more 
involved mean- field potentials [H. Also time-dependent 
simulations are possible using time-dependent density 
functional theory Q. However, hardly any results ex- 
ist, which go beyond mean-field (a rather old example is 
[Toll). Thus this approach works for arbitrary tempera- 
tures if correlation effects due to the Coulomb interaction 
in the system are small. 

Another approach is to start from the von-Neumann 
equation 



-[P-H] 



(4) 



treating the time evolution of the density operator p con- 
taining both the system and the leads. As only the time- 
evolution of the system is of interest, one can trace-out 
the lead degree of freedom in order to obtain the reduced 
dot density operator 



fully reproduced for systems without interactions at ar- 
bitrary temperatures, and the first order generalized mas- 
ter equation schemes, which are recovered in the limits 
of high temperature or large bias. 

In this paper we provide a slightly different presenta- 
tion of our 2vN method in Sec. |lT] in order to highlight 
its structure. Its range of validity and previous results 
are summarized in Sec. IIIIl In Sections IIVIVI we show 
new results demonstrating its applicability to the cases 
of elastic and inelastic cotunneling, respectively. The fail- 
ure of mean-field models, providing a fictitious bistability 
is addressed in Sec. |V]as well. 



II. THE SECOND ORDER VON NEUMANN 
APPROACH 

Using Ho in its diagonal representation (jS]), the tun- 
neling between the states in the leads and the dot reads 
(see Appendix A of 



Pdot(0 = Trioads {p{t)} 



(5) 
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The key task is to find an approximate equation of motion 
for pdot(0 taking into account the coupling Ht- This 
provides (generalized) master equations jTlj . 

To lowest order in Ht electron transitions between 
the quantum dot and the leads are treated indepen- 
dently from each other (here called first order follow- 
ing the number of correlated transitions). Restricting 
to diagonal elements of Pdot{t), one obtains to lowest or- 
der the Pauli master equations 38 1 Furthermore, co- 
herences, represented by nondiagonal matrix elements, 
can be treated in different ways (frequently also denoted 
quantum master equations). The traditional approach is 
the Wangness-Bloch Redfield kinetics [14, JSj which has 
been applied to quantum dots in [l6|, Il7l | While reason- 
able results are typically found, there is a fundamental 
problem due to the occurrence of negative occupations 
as the equations are not of Lindblad type |18l. Lindblad- 
type kinetic equations were derived in ^] , but only 
for the high-bias limit. All first-order approaches entirely 
neglect broadening effects but can treat the interactions 
in the system exactly by a diagonalization of the isolated 
dot Hamiltonian [21| 



The matrix element Ti,a{kat) is the scattering amplitude 
for an electron in the state kal tunneling from the lead 
onto the dot, thereby changing the dot state from the 
state I a) to the state \h). Note that this amplitude van- 
ishes unless the number of electrons in state iV^, 
equals Na + 1 . Here we use the convention that the par- 
ticle number increases with the position in the alphabet 
of the denoting letter. 

A general state vector for the entire system is written 
as \ag) = \a) ® \g), with \g) = \{NMa}) denoting 
the state of both leads, where N^ai S {0, 1}. To 
ensure the anti-commutator rules of the operators we 
use the following notation \g — ka£) = Cf.^(\g) and 
\g + kat) = c\^f\g). The order of indices is opposite to 
the order of the operators. E.g. \g — k'a'i' + kal) — 

- -|g + kat - k'a'l'), 



1.9) = -c^V'^'cLfls) 
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taking into account the anti-commutation rules of the 
operators. To simplify the notation, at is only attached 
to k the first time the index k appears in the equation, 
and in the following it is implicitly assumed to be 
connected with k. Furthermore, we apply the convention 
that X]fc(T(f) iii6ans summing over k and a with a fixed i 
being connected to k in this sum. 



where the states \a) are many-particle states, which may 
be highly correlated. 

In order to treat both broadening and interaction, we 
introduced a second order approach based on the von 
Neumann equation Q which considers both the dynam- 
ics of pdot(i) and higher-order tunneling transitions (2^ . 
This second order von Neumann approach (2vN) bridges 
the gap between the transmission formalism, which is 



The density matrix elements are defined as 



(8) 



where the label n provides the number of electron or hole 
excitations needed to transform g into g' . Examples are 

P^bl-k;ag and P^cl-k-an+k'- We denote the elements as n- 
ehx elements in the following. 
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The von Neumann equation Q gives equations of mo- 
tion of the type 



c.kcri 



E PbL-kTcAk) E PbL+kTb'aik), (9) 

c,kai a,ka£ 



and 



^likTi) as = (p£L,^,,)* and E.P^U.b, = 

T 

rag':bg'—k' 

Performing the sum to the time evolution (jlOp 

some terms do not directly obtain the form $f,J. Here 
we approximate 

[0] 



E ^Nk,lPb'g:bg ~ fk E Pb' g-.bg ~ fk'^b'b^ 



E^cg-fc;c'3-fc = E 'O^^S' ;c'g' ~ (1 - fk)^ 



[0] 



^^-^Pcg-kae-.bg ~ (^c ^ Eb - Ek)p^^g_k.bg 

+ ^Tcb'{k)SN„lP^I;!l.bg-^Pci~k;c'g-kTc'b{k) 
b' c' 

P dg—k~k' \bg 

k'cr'l' b' d 

~ Pig~k-c' g-k'^^c'bik') - P[a-k:a0+k''^baik') 
c' a 

(10) 

Thereby we obtain a hierarchy of n-ehx density matrix 
elements, where the n-ehx elements show a phase rotation 
due to the energy difference involved and are coupled to 
different (n — 1) and {n + l)-ehx elements. In order to 
break the infinite hierarchy we neglect all n-ehx elements 
with n > 3 and obtain a closed set of equations. 

For all density matrix elements we perform a sum over 
all possible lead configurations g and define 
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which are the elements of the reduced density matrix 
Pdot, and 



(12) 



which describe the transitions of electrons between the 
leads and the quantum dot. The particle current from 
the lead ^ into the structure, Jf, equals the rate of change 
of the occupation in the lead. This gives 
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I E ^{E^^we...} 
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I E 3{T4(fc)$W(fc)}. 



(13) 



ka{l),cb 

Thus the $1^^' (fccrf) terms describe current ampli- 
tudes. Note that all 1-ehx terms can be described by 



Similar approximations are done for the p\,g-k-ag 
ments appearing in the equation for the 2-ehx terms. 
This approximation is a factorization of the lead occu- 
pations, which are assumed not to be affected by the 
transition processes. The result is a closed set of dif- 
ferential equations for the reduced density matrix $[,?J,, 
the current elements <i>[,^'(fcCT^) and the similarly defined 

2-ehx terms $[,?b(-/cCT^ + fc'cr'f ;0), (-fccr^- feVT ; 0). 

Defining a discrete set of fc-states, one can set up a col- 
umn vector consisting of all the elements of the density- 
matrix $ — ($M, #[^1, $t^l), where the sub- vectors con- 
tain all the elements of the density-matrix with a specific 
n- value, as well as the complex conjugates of the complex 
elements. The equation of motion for the vector $ can 
be cast on a matrix form 

A ( Eoo Moi 

ih-^ = Mio Ell Mi2 I * = M*. (14) 

V M21 

The submatrices E„„ are diagonal and contain the energy 
differences between the states involved. 

Now we consider the terms in a stationary approx- 
imation, yielding 
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where the iO"*" ensures causality, corresponding to the 
Markov limit for the highest-order elements. 

Inserting the result into Eq. (|14p leads to the matrix 
equation 



ih- 



dt 



$[0] 

Mio Ml J 
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(15) 



where M^^ = -I- M12 (-E22 + iO+) ^ M21 is not di- 
agonal. Due to the fc-dependence it is very laborious to 
express ^^^^ solely in terms of $["1 thereby reducing the 
problem to a generalized master equation. The explicit, 
and completely general, expressions for the equation of 
motion for and $'^1 are given in Eqs. (10,11) of [22| 
and is the main result of that paper. [39] 

The sub-matrices M^g and Mq^ only contain elements 
proportional to the tunneling amplitude T^,a^ while the 
matrix M^^ involves terms proportional to T^^j. Thus the 
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stationary solution of Eq. dTSl) , together with the normal- 
ization ^i'b = Ij ^ill contain all powers of Tba, and 
so will the stationary occupations and coherences. That 
is, the approach does not provide a systematic expansion 
in powers of the tunneling coupling. For the current, 
which is proportional to Ti,a^^^\ all terms up to order 
Tjf^ as well as a class of higher-order terms are taken into 
account. 

Summarizing, the 2vN approach considers the reduced 
density matrix (|lip and the current amplitudes (|12p as 
variables, which are determined by the closed set of dy- 
namical equations (jisp. In deriving these, three approxi- 
mations are applied; (i) only coherent processes involving 
transitions of at most two different fc-states are consid- 
ered, i.e. all n-ehx terms with n > 3 are assumed to 
vanish, (ii) the time-dependence of terms generating 2- 
ehx terms is neglected which corresponds to the Markov 
limit, (iii) the level occupations in the leads, fka-e, is un- 
affected by the couplings to the dot, i.e. the densities in 
the leads and on the dot can be factorized. 



[22| which does not exhibit the oscillations found from a 
time-dependent Green function approach 30]. 

In a recent paper, Jin et al. also consider quantum 
transport in the same spirit as in the 2vN approach by 
keeping correlations between the leads and the dot and 
performing an expansion in the tunneling Hamiltonian 
|3l| . They report a proof that they obtain the 2vN ap- 
proach as a second-order expansion. 



IV. ELASTIC COTUNNELING 

As an example we apply the 2vN method in the elastic 
cotunneling regime for a two-level spinless system. We 
show that the 2vN results agree with a mean-field solu- 
tion and with a scattering result. 

The system is described by the Hamiltonian 



H ^ ^ Eklcl^Ckl + ^ ^Vkincl^dn + h.c 



kin 



(16) 



III. RANGE OF VALIDITY AND COMPARISON 
WITH OTHER METHODS 

For noninteracting systems, Eq. ([3]) gives the correct 
result for arbitrary temperature and bias. In [22j we 
demonstrated analytically that using the 2vN approach, 
Eq. ^ is fully recovered for a single-level system. Nu- 
merically, we also found full agreement for all ranges of 
parameters, including double dots jl^l and the ferrornag- 
netic Anderson model with an applied magnetic field [24| , 
where coherences are of central importance. [401 Thus we 
have strong indications that the 2vN method is able to 
treat transport correctly for noninteracting systems over 
the full temperature and bias ranges. 

For interacting systems and temperatures fcsT ^ F 
or in the high bias limit, the 2vN reproduces the well- 
known results of the methods presented in [20, 25] . This 
demonstrates the correct treatment of charging effects in 
interacting systems. 

For the Anderson model with infinite Coulomb repul- 
sion, the 2vN equations could be solved analytically, see 
[^S*], and the result agree with the diagrammatic real- 
time transport theory in the resonant tunneling approxi- 
mation [2^ , where the onset of Kondo physics is observed 
(see also Eq. (2) of [13]). However, in the Kondo limit 
itself, the model misses the unitary limit and unphysical 
results are found, as the strong correlations between the 
dot and the leads are not properly refiected. 

The validity of the 2vN approach for time-dependent 
problems has not been carefully investigated. As the 
Markov approximation is invoked, it might not be 
valid for strongly time-dependent systems, where non- 
Markovian effects are important due to memory effects, 
which are also relevant when evaluating higher-order mo- 
ments, as e.g., the noise [l^, [2^. As an example, the 
current through a single spinless level was presented in 



with n = 1,2 denoting the two dot states. The first 
term is the Hamiltonian of the leads, and the last two 
terms are the two single- particle states of the dot and 
the interaction between them. The second term is the 
tunnel Hamiltonian with the tunneling amplitudes Vkin- 
Below it is assumed that Vm„ = Xintk, i.e. the 
couplings between both dot states and the lead states 
have a fixed phase factor x^„, and tk is assumed to be 
a real number. The coupling constants are defined as 
TeniE) = 2TTj2^\Vken\^6iE - Ekin) = \xin\^TiE). For 
the 2vN calculations a constant value F for \E\ < 0.95W 
is used, and it is assumed that T{E) = for \E\ > W. 
For 0.95M^ < < W an elliptic interpolation is applied. 
Using the Hamiltonian above, the system can be in four 
different many-particle states, denoted |0), |1), |2), \d) = 
dld\\0), with energies 0, Ei, E2, Ed, respectively. 

The transport is calculated in a setup with Ei <^ //^ 
and E2 + U ^ fie such that the state |1) is the ground 
state. Sequential tunneling processes are blocked due 
to the Coulomb interaction between the electrons, but a 
leakage current due to elastic cotunneling processes can 
occur. Here, the 2vN results are compared with a mean- 
field solution embedded in a nonequilibrium Green func- 
tion framework, and a scattering formalism. 

In the mean-field solution, the interaction term in the 
Hamiltonian is replaced with 



Ud\didld2 



U 



{ [d\di (dld2)+dld2 (d\di^ 



d\di 



■ dldi 



)]}• 



(17) 



where the first [. . .]-bracket is the Hartree term and 
the second the Fock term. The occupations are cal- 
culated self-consistently whereafter the current can be 
evaluated. |4l| 
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For the scattering result, the elastic second-order scat- 
tering rate, 7/^^, is the sum over all processes where an 
electron labelled k' L has been transferred from the left 
contact to the state kR in the right contact, leaving the 
dot state unchanged. This can be realized in two differ- 
ent ways, |1) — > |0) — |1} and |1) — > \d) — > where the 
amplitudes for these two processes are added coherently. 

The rate is calculated as follows ^]: Assume that 
initially the dot is in state |1} and the leads are in a state 
Wli^r)i i-C- the initial state is \i) = \vli^r1) with energy 
Ei = Eiy^^j^ + El. The probability for the leads to be in 
the state vlvris denoted W^^ = W^^^ j as the leads 
are assumed uncorrelated. In the final state an electron 
k'L has been transferred from left to right ending up in 
the state kR, i.e. the final state is \fkk') = c^^Cfe'i|z) 
with energy Ei — E^' l + Eur. According to the T- matrix 
formalism the second-order scattering rates are [3^ 
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Ei — Hq 



■HtV) 



and after some algebra we arrive at 



7fr^ 



r2 



El- E 
X npiE 



E2+U-E 



(18) 



(19) 



where np{E) = [1 + e^/'^s^J^i is the Fermi function, 
and energy-independent coupling constants are assumed 
(the Wide-Band Limit). The rate is found by inter- 
changing L R. Below we only calculate the integral 
for ksT = 0, but note that for finite temperatures the in- 
tegral diverges and a regularization procedure is needed 

For ksT — 0, ^ir = eVbias > and = 0, we obtain 



7ii 



2tt 



dE 



'-LlXRi 



XB.2XL2 



El 



E2 + U -E 



(20) 



and Ji^- — 0. Assuming the state |1) to be almost com- 
pletely occupied (i.e. \Ei\/T < 1 and {E2 + U)/r » 1), 
the second-order elastic cotunneling current is /°i '=otun = 
1 RL 
h^ii ■ 

In Fig.[2]the current versus bias voltage is calculated in 
the elastic cotunneling regime using the 2vN method, the 
mean-field approximation in a Green function framework, 
and the second-order scattering method. For the latter, 
the calculation is for vanishing temperature, while the 
other calculations are for ksT = F/IO. For Ei = — SF 
almost perfect agreement between all three methods is 
found, while deviations between the scattering method 
and the others occur for Ei = — 2r, which is most likely 
due to the fact that the state |1) is not fully occupied in 
this case as assumed in the derivation of the second-order 
scattering rate ([TO)) . 



FIG. 2: (Color online) The /V^-characteristics for two different 
values of -Ei calculated using the 2vN approach, mean-field 
(MF) nonequilibrium Green functions, and, finally, the scat- 
tering cotunneling current 7=1 1^"'"" = i(7ii'^ -711^) , where 
the latter is calculated for ksT — 0. The chemical potentials 
are fj.L = eVbias and fiR = 0, and the phase factors are all 
equal xin = l/\/2, i-e. Fin = r/2. The other parameters are: 

E2 = ler, u = 2or, fcsT = r/io and w = 5or. 



In summary, from a comparison with both a mean- 
field solution and a scattering formalism, we have shown 
that the 2vN method is able to quantitatively describe 
elastic cotunneling processes even for temperatures much 
lower than the energy scale set by the coupling to leads, 
/cbT<F. 



V. INELASTIC COTUNNELING 

For higher bias voltages the agreement between the 
2vN approach and the mean-field solution is assumed to 
be less perfect, as the latter can lead to bistable solutions 
[3^ , which are not present in a generalized master equa- 
tion approach. Fig. |3] shows a comparison between the 
2vN approach and the mean-field solution over the full 
bias range for Ei = — 2F, with the rest of the parameters 
as in Fig. O 

Considering first the 2vN result, the curve shows in- 
creased current when the bias matches the energy differ- 
ence between the levels, eVbias = I 1 ~ ^2 1 = 18F. This is 
due to the onset of inelastic cotunneling [33] , which leads 
to a population of the excited state, |2). After the inelas- 
tic cotunneling process, additional cotunneling-assisted 
sequential tunneling processes can occur [H,!!!]. Finally, 
at eVbias = E2 + U = 36r sequential tunneling through 
the upper level becomes possible. The value for the cur- 
rent is consistent with the master equation result (3F/8?i) 
plus additional cotunneling through the lower level. 
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FIG. 3: (Color online) The /V^-characteristic over the full bias 
range for i5i = — 2r calculated within the mean- field approx- 
imation and using the 2vN approach. The other parameters 
are as in Fig. [2] The arrows indicate the direction of the bias 
sweep in the mean-field calculation. 



In contrast, the mean- field solution misses the onset of 
inelastic cotunneling and instead a bistable behaviour is 
observed, and it also overestimates the current after the 
onset of sequential tunneling. Both is due to an insuf- 
ficient description of the Coulomb interaction and em- 
phasizes the need for a method which can describe both 
higher-order tunneling processes (as elastic and inelastic 
cotunneling) and many-particle interactions. 



VI. SUMMARY 

The second order von Neumann approach provides 
a quantitative description of transport through nanos- 
tructures for arbitrary bias and temperatures above the 
Kondo temperature. In the low-bias regime, the elastic 
cotunneling current is in very good agreement with both 
the second-order scattering result and the mean-field so- 
lution, even for temperatures much lower than the energy 
scale set by the coupling to leads. Inelastic cotunneling is 
also well captured by the 2vN approach, whereas a mean- 
field solution within a nonequilibrium Green functions 
framework shows an artificial bistability in this regime. 
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